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Heavy ion collisions at RHIC are well described by the (nearly ideal) hydrodynamics for average 
events. In the present paper we study initial state fluctuations appearing on event- by- event basis, 
and the propagation of perturbations induced by them. We found that (i) fluctuations of several 
lowest harmonics have comparable magnitudes, (ii) that at least all odd harmonics are correlated in 
phase, (iii) thus indicating the local nature of fluctuations. We argue that such local perturbation 
should be the source of the "Tiny Bang", a pulse of sound propagating from it. We identify 
its two fundamental scales as (i) the "sound horizon" (analogous to the absolute ruler in cosmic 
microwave background and galaxy distribution) and (ii) the "viscous horizon" , separating damped 
and undamped harmonics. We then qualitatively describe how one can determine them from the 
data, and thus determine two fundamental parameters of the matter, the (average) speed of sound 
and viscosity. 

PACS numbers: 



I. INTRODUCTION 

Starting the introduction, let us note that the issues to 
be discussed in this paper are somewhat similar in nature 
to current trends in cosmology of the last decade. While 
the very existence of cosmic microwave background radi- 
ation had dramatically confirmed the existence of the Big 
Bang already in the 1960's, it is its more recent observa- 
tions which made cosmology a really quantitative science. 
Small temperature fluctuations on top of the overall Hub- 
ble expansion have been seen on the sky. Their angular 
size and the magnitude of various harmonics tell us, in 
surprising detail, about sounds propagating the Universe 
at the plasma neutralization time, providing precise tim- 
ing of the cosmological expansion. 

Experimental data obtained in heavy ion collisions at 
the Relativistic Heavy Ion Collider (RHIC) has found 
the "Little Bang", a hydrodynamical explosion driven 
mostly by pressure of the new form of matter, Quark- 
Gluon Plasma. Their experimental data for radial and 
elliptic flows has been compiled in the so called "white pa- 
pers" of all the experiments in 2004, and compared with 
predictions of relativistic hydrodynamics. Very recent re- 
sults from the Large Hadron Collider (LHC) on elliptic 
flow [1 also turned out to be in agreement with hydro- 
dynamical predictions, suggesting that QGP remains a 
good liquid even at LHC (see e.g. [2J). (For clarity, 
we mean hydrodynamics complemented by the hadronic 
cascade for a correct account of the hadronic stage, see 
[3]-[5]. Models which do not do that and use rather ar- 
bitrary freezeout don't get the energy dependence right, 
neither for RHIC nor for LHC.) Dissipative effects from 
the QGP viscosity provide only small corrections at the 
few percent level, see [6-8 . So, by now, we have a good 
quantitative description of the "Little Bang" . 

This paper is however not about it, but about phenom- 
ena which we will call collectively "The Tiny Bang" . In 



the hydrodynamical studies mentioned above one ignored 
any fluctuations in the system, using average smooth 
initial conditions, possessing the symmetries of the av- 
erage collisions. (For example, central collisions are 
azimuthally symmetric, with vanishing nonzero Fourier 
components of the flow. Non-central collisions were as- 
sumed to produce reflection symmetric x — » — x distribu- 
tions at midrapidity, thus only even harmonics are to be 
nonzero.) Yet local perturbation of those should produce 
extra excitation of the "Little Bang" , that should not in 
general respect such symmetries. 

Such perturbations of the average explosion can come 
from at least two different sources. The one which we 
will study in this paper is due to quantum fluctuations 
in the wave function of the colliding nuclei, which creates 
"bumpy" distributions of matter, for any collision, which 
one can decompose into a smooth average one plus local 
perturbations. 

Another one, to be studied in subsequent papers of this 
series, are created by the energy deposited by jets propa- 
gating through the medium. It has been recently dramat- 
ically shown by ATLAS collaboration [9] that even jets 
with energy above 100 GeV deposit large part of their 
energy, and sometimes all of it, into the medium. The 
first ideas were to look at the resulting perturbations in 
the form of a Mach cone [10] , driven by the view that the 
energy is deposited more or less homogenously along the 
jet path. However more recent developments of the the- 
ory, based on AdS/CFT, have lead to the view that the 
deposited energy grows as the cube of the distance trav- 
elled by the jet, AE ~ L 3 , with significant deposition 
near the endpoint. Thus one may think of the second 
kind of the "Tiny Bangs", this time occurring roughly 
half time between the beginning and the end of the "Lit- 
tle Bang" . 

The smallness of the perturbation amplitude, with re- 
spect to the local density of ambient matter, would sug- 
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gest the appearance of divergent sound waves, see FigjT] 
Similar to the circles from a stone thrown into a pond, 
hydrodynamics tells us that initial perturbations should 
become moving waves, with basically nothing left at the 
original location at later time. This is the picture we are 
going to work on in this paper. 




FIG. 1: A sketch of the transverse plane of the colliding sys- 
tem: the two concentric circles are the nuclear radius (inner) 
and the final radius of the fireball (outer). The black spot in 
the upper figure is an extra density at formation time, due 
to initial density fluctuations in the collision. The small per- 
turbation becomes a circle of a sound wave, which will be 
stopped at freezeout at thick circle (distorted by radial flow). 
The part outside of the fireball does not exist: the correspond- 
ing matter will actually be placed near the edge of the fireball 
(thick dashed line). The whole perturbation is enclosed in a 
sector between the two thin lines with arrows. 



An alternative idea, of randomly fluctuating shapes of 
the produced initial fireballs, has resulted in an approach 
in which different angular harmonics of that distribution 
are treated separately. The realization that even central 
6 = collisions may have some fluctuating ellipticity has 
lead to the discussion of the elliptic flow event-by-event 
fluctuations, see [IT] and many subsequent works. The so 
called "triangular flow" related to the 3-ed harmonic of 
the flow has been recently studied by Alver and Roland 
[12], with several groups working in this direction now. 

The main difference between our approach and that is 
that we treat such fluctuations not as independent noise 
in different harmonics but as certain local perturbations, 
resulting in certain evolving sound fronts, reaching cer- 
tain size, shape and diffusivity by the moment of freeze- 
out. On one hand, one may argue that as soon as all the 
perturbations are small and the equations are linear, it 
is not important if one expands in harmonics before or 
after the solution of hydro equations. And yet we believe 
that our approach is not only more intuitive and pro- 
vides better insights, but it corresponds to the physical 
nature of the fluctuations in question, as higher angular 
harmonics are in fact correlated with each other as we 
show below. Furthermore, we think that our approach is 



in more direct connection with the fundamental scales of 
the problem, to be detailed below. 

Having outlined the main ideas of this paper, let us 
discuss a bit more some recent relevant papers. In the last 
few years RHIC experiments have focused more on two 
and three-particle correlations, which revealed a rather 
rich phenomenology of correlations. We are not however 
able to review those, and we will refer to data as we 
develop the theory more. 




FIG. 2: Extra particle distribution in azimuthal angle relative 
to position of the initial perturbation. Fig. (a) from [13] has 
six curves, from the most narrow to wider ones, correspond 
to the radius of the circle 1,2,3,4,5,6 fm, respectively: the 
last (blue) with a maximum away from the original position, 
corresponds to the sound horizon. The original spot position 
is selected to be at the edge of the nuclei and the observed 
particle is at pt = 1 GeV. 

Fig.(b) from [14] . for a particle with p t > 2 GeV 



The propagation of sound on top of the fireball has 
been discussed by J.Casalderrey-Solana and one of us in 
[T6] . In that paper the fireball expansion was modelled 
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by the Big-Bang-like overall expansion of the space, with 
the same Friedman-Robertson- Walker metric as used for 
cosmology. The focus of that paper was the effect of time- 
dependent sound velocity, especially if the phase transi- 
tion is 1st order and it can vanish at some interval of T. 
The interesting finding was a creation of the secondary 
- and convergent - sound waves. This idea was further 
discussed in [13] in connection with the "soft ridge" is- 
sue, but with the conclusion that if the current lattice 
data on the speed of sound is correct, the effect of the 
reflected wave is too small. 

In the same paper [13 j it has been found that the usual 
(unreflected) sound propagation should produce charac- 
teristic "two-peak events", with the angle between the 
peaks reflecting the sound horizon and numerically being 
about 1 rad, see the blue curve in Fig. [5|a). Andrade et 
al [14] have independently come to the same conclusion, 
see Fig. j2^b) taken from their paper. Note that their 
peaks are at the same angle: they are more pronounced 
simply because the particle p t for which it is plotted is 
higher. Note also a dip near zero, indicating that nothing 
remains at the original location of the bump. Andrade 
et al have further pointed out that the two-peak events 
lead to a three-peak correlation function shown in Fig. 
[3] with the side peaks now at twice larger angle ±2 rad, 
now on the "away" side from the peak. (Note that the 
central peak is about twice larger than the peripheral 
ones: this is because it corresponds to events 11 and 22 
while the other peaks are 12 and 21 combinations, from 
the peaks 1 and 2 in the two-peak distribution. ) 

This observation explained what has been found ear- 
lier, in the "event-by-event" hydrodynamical studies by 
the Brazilian group (see [18] and references therein). The 
Brazilian group has used the initial condition from string- 
based model developed by Werner and collaborators. Re- 
cently this group also studied event-by-event hydrody- 
namics [15] and claimed that their correlation function 
describes the data on two-particle correlations. They 
trace their origin to multiple "bumps" in the initial distri- 
bution that lead to the development of what they called 
"fingers". There are about O(10) of them per event, each 
corresponding to rather narrow angles. Presumably the 
"event-by-event" complicated pictures are more or less 
linear superpositions of the independent two-peak struc- 
tures from each bump discussed above, because the sound 
waves (as any Goldstone modes) hardly interact with 
each other. We will discuss all those issues in detail: and 
for now let us only note that all these works have not 
yet included the key phenomenon to be studied, namely 
viscosity. It is easy to estimate that even the minimal 
conjecture viscosity n/s = 1/4tt would significantly mod- 
ify the amplitudes of those fluctuations. 

The most striking outcome of all these works p~2] [15] 
[T8] is that they have obtained the "Mach-cone-like" cor- 
relation functions, with the away-side ridges at ±2rad 
(this in the other hemisphere from the trigger) without 
any jets involved! This created large enthusiasm, and 
attempts to disregard hard collisions and jets, even for 



the events with trigger hadrons of few GeV p t . How- 
ever, there is no doubt that such particles would not be 
produced from hydrodynamics and are still due to hard 
collisions producing jets. 

Although this work is in many respects a continua- 
tion of the paper [13 , we will focus on sounds and post- 
pone the development of another important idea pro- 
posed there. In spite of its phenomenological success, 
there are good reasons to think that the hydrodynami- 
cal description of the RHIC evolution may be incomplete 
because in the so called M-phase of the collision elec- 
tric gauge fields should remain unscreened for the time 
of evolution. It has thus been proposed to upgrade hy- 
drodynamics to (dual)magnetohydrodynamics, including 
these unscreened electric fields in the stress tensor. This 
changes the nature of the propagating modes, and it also 
makes possible the survival of the "flux tubes" in the 
"QGP corona" . We hope to return to those effects later. 




FIG. 3: (from [14]) the two-particle correlation function from 
the two-peak events, for particles with transverse momenta 
indicated on the figure. 



II. SETTING UP THE PROBLEM 
A. The main scales of the problem 

Before going to specifics, let us formulate the problem 
in a more general form, which in a way connects it to the 
Big Bang fluctuations. 

Two generic scales are (i) the macroscopic scale R and 
the microscopic scale Z, being in the relation 

1<^R (2.1) 

which ensures such macroscopic tools as thermo and hy- 
drodynamics to work. 

The macroscopic scale R is the size of the fireball in 
heavy ion collisions and the curvature scale a in the Big 
bang. Note that both are in principle time dependent, 
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demonstrating the expansion of the system. However, 
while a(t) changes by many orders of magnitude, the fire- 
ball size increases rather modestly, e.g. from 6 to 8-9 fm 
at its maximal size, for AuAu collisions at full RHIC en- 
ergy. 

The microscopic scale / is the mean free path for weakly 
coupled systems (weakly coupled QGP or hadronic gas): 
in the case of strongly coupled QGP (sQGP) it is just 
the inverse temperature I = 1/T. For AuAu collisions at 
full RHIC energy I changes from .5 to about 1 fm, from 
initial to hadronization time. Thus the large parameter 
R/l = O(10) in the region we use hydrodynamics. 

Now let us define two new scales. The first is the sound 
horizon 



r f 



drc s (r) 



(2.2) 



where the integral is taken from the formation to freeze- 
out time. At the freezeout the waves just stop where they 
are, and the matter is split into independent particles. 

It is the same idea as suggested by Sunyaev and Zel- 
dovich for the Early Universe [T7] : the initial perturba- 
tion (say higher density at some point) creates a sphere 
of such radius, at which the density is a bit higher than 
the average: when galaxies are formed they are corre- 
lated with that sphere and thus is observed today in their 
correlation function. It turns out that in the Big Bang 
this produces a "standard ruler" , which is today of about 
H s w 150 Mps, observed in galaxy's distribution and in 
CMB correlations. 

One of the main issues discussed in this paper is 
whether any manifestation of the sound horizon scale can 
be observed in the Little Bang. For example, one may 
think of angular correlations with angles 



A0: 



2H S 
R 



(2.3) 



or angular harmonics with m ~ 1/A0. In the cosmology 
such angular momentum is I ~ 200. Going ahead of 
ourselves, we will show that in our problem of the Little 
bang, for AuAu collisions at RHIC, we will deal with 
m ~ 3. 

The second scale is not important in cosmology and 
we would like to call it "the viscous horizon scale 11 R v . 
Its verbal definition is that it separates the wavelengths 
of the sound which are and are not dissipated by the 
viscosity effects. The smooth fireball and fluctuations 
are described by 



6T n 



(2.4) 



The textbook dispersion law for the sound, including the 
viscosity term, is 

uj = c s k -— (2.5) 

2 3s T v J 

After that Fourier transform puts it into momentum 
form, after which one can solve the time dependence us- 
ing the momentum-dependent dispersion relation as well 
as the imaginary part induced by viscosity. 






FIG. 4: Two upper picture correspond to initial time t = 
0: the system has almond shape and contains perturbations 
(black spots). Two lower pictures show schematically location 
and diffuseness of the sound fronts at the freezeout time £/. 
The arrows indicate the angular direction of the maxima in 
the angular distributions, 2 and 3 respectively. 



(One may add bulk viscosity to this expression as well, 
but we keep the shear viscosity for now, assuming it is 
dominant.) 



ST^i) = exp 



The spectrum of the original t — perturbations have 
harmonics of the so called "saturation scale" Q S1 which 
is for RHIC of the order of Q s ~ lGeV. Even if one 
takes the minimal viscosity n/s = 1/47T, by freezeout 
t ~ 10 fm/c this exponent gets very large, damping such 
fluctuations to un observably small magnitude. Only the 
harmonics below the new viscous survival scale k < k v 
would survive, which is determined from the condition 
that the exponent above is less than 1 



k v — 



2tt 
R v 



l3Ts 

27/7? 



200MeV 



(2.7) 



(the number comes from an estimate rj/s = l/4rr,T ~ 
200 Me V, 7/ ~ 10 fm/c). 

What it means is that sound perturbation spheres 
(cylinders, cones etc) would not have the width of the 
original fluctuations, but they get significantly widened, 
with the width of the order of l/k v . Note that while the 
radius of the spheres increases linearly with time ~ t, this 
width increases only as t 1//2 , which means that although 
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the spheres become more diffuse, they are also relatively 
sharper and sharper as time goes by. 

Let us finish this section by pointing out the main aims 
of this investigation. By observing propagating sound 
perturbations one would like to measure the two scales, 
H s and l/k v , experimentally, defining two key hydro- 
dynamical parameters - the speed of sound and viscosity. 
The way to do so is to change the geometry of the collision 
(by centrality) and the size of the nucleus (by changing 
the beam A), and by observing all the harmonics of the 
flow. The amplitudes of the higher harmonics, dampened 
by viscosity, if measured, would provide an independent 
measure of the viscosity. 

For central AuAu collisions at RHIC the hierarchy re- 
lation between all those four scales is 

R> H s > R V >1 (2.8) 

As some representative numbers let us mention 
8,4, 2,0.3 /m, respectively. The observation angle of the 
"peaks" is 

A0«2^ (2.9) 

Their angular width is 

5<j> = R v /R (2.10) 

so harmonics larger than m > R/R v can hardly be found, 
as they are dissipated. 

However for mid-central collisions the width (short 
size) of the "almond" R x becomes comparable with H s , 
and for more peripheral ones R x < H s . As a result, 
one expects the sound to traverse the whole fireball and 
deposit some amount of (entropy) density to its side op- 
posite to the original fluctuation, as it speeded up with 
the flow. In this case one expects 3-peak events: see Figj4] 
for explanation. 

So, one of our suggestions to experimentalists is to lo- 
cate the centrality at which such transition occurs. 

III. THE INITIAL STATE FLUCTUATIONS 

A. Generalities 

Let us start with a comment on what we would call 
the "initial state" . This term is currently used in at least 
three different settings: 

(i) The wave function of the colliding nuclei, expressed 
either in terms on the nucleons (their positions in the 
transverse plane just prior to the collisions) or in terms 
of partonic degrees of freedom (positions and longitudi- 
nal momenta). Another version of it is the "Color Glass 
Condensate" (CGC) described as an ensemble of classi- 
cal gauge fields. 

(ii) The state just after the (Lorentz contracted) nuclei 
passed each other. It is either the partonic state, includ- 
ing partons newly produced in a collision; or the so called 



GLASMA, in the classical field description. 

(hi) The state after approximate equilibration is reached, 

so that macroscopic (hydrodynamical) description can be 

started. 

It is the last one which we mean in this work, as we 
would apply hydrodynamics as a tool, translating prop- 
erties of the initial conditions into the final state observed 
in the experiment. Therefore our "initial state" should 
correspond to about one unit of the relaxation time af- 
ter the actual collision, or numerically at a proper time 
of the order of 1/2 fm/c. Thus the inhomogeneity of the 
initial wave functions should be already smoother than at 
time zero, by this (so far poorly understood) relaxation 
process. 

As we detail below, this state will be described by some 
"average" or zeroth-order shape of the fireball (depending 
of course on the impact parameter, the colliding nucleus 
and the collision energy), plus "fluctuations" character- 
ized in the first order by an ensemble of small pertur- 
bations of the average shape described by Fourier coef- 
ficients and phases {e n ,^ n }. Generic expressions would 
include the zeroth order ensemble- average deformations 

< e n > and deviations which have no average but fluctu- 
ations 5e^ =< e 2 n > — < e n > 2 . 

The simplest situation, happening for the second har- 
monics and sufficiently peripheral collisions, is that the 
average is much larger than the fluctuations, < >^> 
5e2. If so, one may assume Gaussian form of the fluctu- 
ations with the width given by Se2. But the situation 
is quite different for near-central collisions, for which 
both terms in e n come from fluctuations. Their dis- 
tribution are obviously non- Gaussian because they are 
all positive by construction. So one must introduce and 
study their higher powers or cross-correlations (such as 

< Se ni Se n2 Se n3 > and their generalizations with certain 
combination of phases (see below). All of those should 
in principle be provided by the "initial state models" , of 
which we select Glauber model as the simplest example. 

The separation of the initial state fluctuations from all 
other fluctuations (e.g. fluctuations during the hydro- 
dynamical evolution, hadronization and the freezeout) is 
possible because of the fundamentally different number of 
relevant degrees of freedom defining their magnitude. As 
we will detail in the next section, the so called Glauber 
fluctuations due to various number of "wounded" (or par- 
ticipant or interacting) nucleons are of the order of 




where the number of the participant nucleons N p ~ 
0(100), being limited from above by the total nucleon 
number 2 A ~ 400. 

Further fluctuations are determined similarly, but with 
the number of participants N p substituted by the much 
larger number of partons involved, or the total multiplic- 
ity Nh a drons °° 10 4 (for RHIC and LHC it is factor 2 
different). That is why one may, to certain accuracy, 
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ignore all later-time fluctuations and assume that ob- 
servable fluctuations in particle spectra and correlation 
functions are one-to-one translated from the initial state 
ensemble. Thus we use hydrodynamical equations as a 
fully deterministic tool, by itself producing no random 
numbers at all. 

Furthermore, for near-central collisions all Se n are 
small, of the order of several percents. So, independently 
of their possibly complicated distributions and cross cor- 
relations, the hydrodynamics applied in linear approxi- 
mation should be quite reliable tool. Thus hydro equa- 
tions can be linearized and the linear response coefficients 
Sv n /Se n calculated. If so, it does not matter what the 
actual magnitude of the deformation Se n is. Also the 
linearized perturbations do not interact with each other. 

Although we will focus on those calculations in our 
next paper, let us note here two things. One simple fact is 
that while angles ip n of the fireball deformations indicate 
the maxima of the distribution (the corners of triangle, 
square and other polygons), hydro flow goes along their 
sides. Therefore the observed flow angles £ n are rotated 
from the deformation angle as follows 



(3.2) 



Our second comment is that higher harmonics n are sup- 
posed to become oscillatory in time, displaying acoustic 
sound properties. Naively one may think that at freeze- 
out this leads to their random phases, and thus those can 
be ignored. This is however not true, as only the coher- 
ent sum of all harmonics with their correct phases will 
reproduce a propagating sound wave. Finally, if needed, 
nothing prevents one from solving full nonlinear hydro- 
dynamics (without linearization). In fact it is done by a 
few groups devoted to "event-to-event hydrodynamics". 
If done one can also calculate the nonlinear effects, e.g. 
appearance of the 4-th harmonics in spectra coming 
not from 5 64 but from ((^2) 2 . 



B. Fluctuations in the Glauber model: the 
amplitudes 

Our "Glauber model" is a bit different from that used 
widely by experimentalists. Both assume that initial 
state fluctuations originate from the nuclear wave func- 
tions. The "usual Glauber" uses randomly placed coordi- 
nates of the individual nucleons in the nuclear wave func- 
tion. However, the nucleons themselves are complicated 
objects and their interactions are also strongly fluctuat- 
ing: since there are studies of that we decided to include 
this source of fluctuations also. This changes numbers 
a bit, but was found not to be important for any of the 
qualitative conclusions to be reached. 

The nucleon fluctuations we included via the fluctu- 
ating NN cross sections which are to a certain degree 
known and studied via diffraction, see [19] for the de- 
tails and earlier references. Naively, from the well known 



fact of a nucleon being made of quite a large number of 
partons one might conclude that those fluctuations are 
small, 0(l/N partons ): but this is not the case. In our 
simulation we have assumed the cross section o~nn to be 
the random Gaussian variable with the variance 



wnn — 



< 0~NN > 2 



< 0~NN > 2 



i 0.25 (3.3) 



First, like in [12], we simulate a large ensemble of col- 
lisions and calculate the magnitude of the e n for several 
lowest harmonics (up to 6). Their definition is via the 
Fourier expansion for a single particle distribution 



fW = 7^ (l + 2^e n cos( 

\ n 



n(^-V„)) (3.4) 



where the e n are the participant anisotropics and the ip n 
are the angles between the x axis and the mayor axis of 
the participant distribution. 

The participant anisotropics are calculated from 



\J (r n cos (ncj))) 2 + (r n cos {n(j))Y 
(r n ) 



(3.5) 



This expression is calculated in the center of mass of the 
participant nucleons for each event. Therefore the dipole 
moment n = 1 made out of the average coordinates 

< x >= (r cos ((/>)) =0, <y >= (r sin ((/>)) = (3.6) 

are zero by definition. 

The 2-d shape of the event can in principle be ex- 
panded in the double Taylor series in x,y or in double 
series over moments r™cos{n<ft),r m sin{ri(j)) with integer 
m, n. An even better definition would be to follow the 
customary statistical trick and write the distribution as 
the exponent containing the "generating function" of the 
angular dependence expended in harmonics 

P = F 1 (r)exp(F 2 ), F 2 = ^ r n e n cos{n{cf> - (3.7) 

n>0 

In this way the positivity of the distribution function, as 
well as inclusion of trivial higher order effects are ensured. 

Since the dipole m = n = 1 term is zero by construc- 
tion, we define the first odd deformation e\^\ using the 
term of the expansion m = 3, n = 1 which appears to- 
gether with the triangular deformation m = n = 3 



yj (r 3 cos (</>)) 2 + (r 3 cos ((j))Y 



(3.E 



The anisotropics calculated in this way are plotted in 
figure [5] for n = 1,6. The plot shows that the eccen- 
tricity has the largest value for the well known elliptic 
deformation 62 and a nonzero value of triangularity 63, 
in agreement with the results reported in [12] . Note that 
for the near-central collisions N part > 300 the elliptic 
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FIG. 5: Average anisotropies (upper plot) and their variations 
(lower), as a function of centrality expressed via the number 
of participants N par t 



deformation is no longer dominant, and it is also due 
to fluctuations. This conclusion becomes evident as one 
looks at the lower plot in Fig|5j which shows the varia- 
tions of these e n . 

One observation coming from these results is that all 
other deformations (except for ei, small because the "true 
dipole" remains zero) are all comparable, ranging from 
0(1/10) for central collisions to 0.3 -0.5 for most periph- 
eral ones. While in the central bins these perturbations 
can be considered small and treated as Gaussian random 
variables, it is clear that for most peripheral bins (when 
the number of participants is smaller) the fluctuations 
are large and thus must be non-Gaussian. 

Another consequence is that there is absolutely no 
ground to single out 63: in fact both 64 and 65 are larger 
that 63 and cq is about of the same order as 63. 

The last point is that their variations (the lower plot) 
are all comparable to the magnitude. Yet the definition 
of deformations are such that they are always positive, for 
each event. This is one more reason that the amplitudes 
cannot have Gaussian distribution, deviating from it at 
least for the smallest values. 
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FIG. 6: Distribution of the angles ip n for the first three har- 
monics, the centrality bin used is 100 < N par t < 300 



C. Fluctuations in the Glauber model: the angles 

The angles ip n are defined by: 

(r n sin (n<j))) 



tan (nip n ) = 
and to calculate ipi we use: 
tan (ipi) — 



(r n cos (n(j))) 



(r 3 sin (0)) 
(r 3 cos {(j))) 



(3.9) 



(3.10) 



Using these expressions we obtain the distribution of the 
^ n 's for the first six harmonics as shown in Figs[6j[7| In 
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FIG. 7: Distribution of0n for the harmonics 4-6, same cen- 
trality 



order to better understand the behavior of these angles 
we will now study their correlation. 

(Note that our angle definition is different from the one 
by Alver- Roland [12]: we do not include extra phase tt/u 
between the flow and deformation directions, see below.) 

Let us comment on these distributions, starting from 
the even ones. 

The most obvious one is a distribution of the sec- 
ond (elliptic) harmonic: as seen in Fig{6] the angle tp2 
is strongly peaked at 7r/2, corresponding to an elonga- 
tion of the system in the y direction, as of course one 
expects from the overlap "almond" of two nuclei. The 
distribution of the 4-th angle ip^ in FigjT] shows peaks 
at angles and tt/2: but since quartic symmetry of the 



4-th harmonics it simply means that the maxima of the 
distribution tend to be aligned with the coordinate axes 
x and y. The distribution of the 6-th harmonics is dif- 
ferent: it is peaked at the angle 7r/6. This means that it 
has no maximum at x direction but rather in y. In con- 
clusion, all even harmonics are strongly correlated with 
the reaction plane, all of them producing maxima along 
the y (out-of-plane) direction. 

The distribution of the angle ipi is nonzero at all angles, 
which means it is not very strongly correlated with the 
reaction plane. It has two maxima, at ±y directions, to 
be called "tip" fluctuations. Although the contribution 
from angles 0,7r or x-directions is about twice smaller, 
it also makes an important contribution: we will call it 
"waist" fluctuations. Note that while the area of the 
"waist" is larger than "tips": and yet its contribution is 
smaller. 

The distribution over ^3, ^5 in these figures looks com- 
pletely uncorrelated with the reaction plane. (This fact 
has also been noticed in [12 1 and by others.) However, 
further scrutiny shows that they are in fact well corre- 
lated with ipi, see Fig. [8] (in which we included points 
repeated by periodicity). The distribution can be crudely 
characterized by some "bumps" plus "stripes" connecting 
them. 

The interpretation of the "bumps" is that all of them 
correspond to events with additional "hot spot" at the 
"tips" of the almond. It is a very natural place for maxi- 
mal fluctuations, for two reasons. First, this is where the 
participant density in both nuclei is near zero. Second, 
because of the factor r 3 they have larger weight than all 
other places. 

There are two kinds of "stripes" , with positive and neg- 
ative slope in Fig. [8] The latter one simply follow from 
distribution, while the former one is indeed a nontriv- 
ial correlation between the angles whose origin we cannot 
explain. We will continue to discuss its manifestation a 
bit later. The correlation of ^5 with ipi is very similar. 
The "bumps" at ^5 — ipi ~ again mean ±y the direction 
or the "tips" . The plot has similar "stripes" . 

Going a bit ahead of ourselves, let us study the "reso- 
nant" combinations of angles, as well as angles and am- 
plitudes. As we explain below, those particular combi- 
nations of the amplitudes and phases or two harmonics 
are 



2 C08(ni^ni -^2^n 2 ) > 



(3.11) 



especially in the case when ni,n2 differ by two units. 
We have studied two first examples of the kind, with odd 
harmonics 1,3,5. 

One interesting distribution, shown in Figj9j is that 
over the cos term itself, for the particular combination of 
the 1-3 phases. It consists of two clearly different parts: 
a very narrow peak near —1 and wide fiat distributions 
between -1 and 1. This plot demonstrates qualitative 
feature of the phase distribution which was pointed out 
above. One explanation of the peak near -1 (the an- 
gle combination is tt) comes from the fluctuations at the 
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FIG. 8: Scatter plot of the 1^3 vs 1^3 — ipi (above), and of the 
ips vs ip5 — ipi (below) , the same centrality 
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FIG. 9: Scatter plot of cos (3^3 - V>0 
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"tips" of the almond, when both tpi and tps are strongly 
correlated close to tt/2. However, the second interesting 
correlation seen as "positive slope lines" in FigjS^i be- 
cause for them ipi — 3^3 = tt as well. A similar situation 
happens for other odd harmonics. 



The average value of the combinations (3.11) for 1-3 
and 3-5 harmonics as a function of centrality are shown in 
FigjTol All values are negative, as the sign is dominated 
by a peak in cos near -1: the other component more or 
less averages out. We thus conclude that experimental 
measurements of the amplitude of such correlations, with 
the magnitude and the sign, will be especially sensitive 
to the "almond tip" fluctuations. 

Summarizing the observed pattern: we have found that 
all odd angles are well correlated with each other, form- 
ing the "stripes" and "bumps" shown in Fig|8] The fluc- 
tuations and correlations seem to be stronger from the 
"tips" of the almond. 



-0.04 



-0.06 



-0.08 



FIG. 10: Correlators (tit^cos (ipi — 3^3)) (top) and 
(6365 cos (3^3 — 5^5 )) (bottom) as a function of the number 
of participants. The error bars are omitted since they are 
smaller than the dots. 



D. Comments on other initial state fluctuations 

So far the only source of fluctuations included has been 
(i) the coordinate part of the nuclear wave function pre- 
scribing the nucleon positions in the transverse plane, 
and (ii) the event-by-event fluctuation of the NN cross 
section. We found that the former effect dominates and 
the latter only provides small corrections. 

While other sources of fluctuations clearly are subject 
for future studies, we still provide some comments on 
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those. 

One important type of "initial state" fluctuations is of 
course hard part on scattering events, resulting in jet pro- 
duction. The rate of those very strongly depend on the 
exact definition of the cutoff beyond which the momen- 
tum transfer involved is characterized as "hard". There 
are vastly different opinions on where this boundary is 
theoretically, and experimentally it depends on whether 
such events are triggered by single large-p^ hadron or by 
some jet-finding algorithms. Jet production and quench- 
ing is of course of high interest, but those should be stud- 
ied only in small fraction of all collisions selected by sep- 
arate triggers. For global fluctuations those can safely 
be ignored, as the probability of "hard" events is smaller 
than that of the fluctuations we study. 

In the Glauber approach that we used, the local den- 
sity of produced matter is assumed to be simply pro- 
portional to the local density of all participant nucleons, 
Np(Ai) + N P (A2). However, when this density is high 
enough, it has been argued that the so called "satura- 
tion" phenomenon should take place, because of parton 
absorbtion processes in the wave function. Other ex- 
pressions for local matter density have been proposed, 
e.g. ~ min(N p )ln(max(N p ) /min(Np)) by Kharzeev, 
Nardi, Levin, where min and max refer to the smaller and 
the larger of the two. Those are typically the amplitudes 
of the fluctuations. 

A well known approach to their description is the so 
called "glasma", which calculates those color fields from 
random color charges of the leading (larger x) partons of 
the two colliding nuclei. Asymptotically (in a very large 
nuclei or at very high energy) McLerran and Venugopalan 
[24] argued that at a particular location in the transverse 
plane the color charges of partons must be uncorrelated 
because they come from different nucleons. Therefore it 
is usually assumed that their color charges fluctuate as 
random variables. If so, the resulting fluctuations are 
small, as the total number of partons is very large. 

Application of such ideas usually keep the average val- 
ues of those, such as A 1 / 3 or so. The point of our 
comment is a warning, that such simplified ideas can- 
not be used for determining the fluctuations. It has been 
known for a long time, that while at very small x the 
partons, mostly gluons, become numerous, they are still 
very tightly correlated in the transverse plane. The size 
of the "gluonic spot" in the nucleon has been known for 
long time from diffractive form factors, and in more re- 
cent form, from HERA photon diffraction into J/ip. This 
spot is small, and therefore bright. As documented e.g. 
in [21] (their Fig. 23), the gluon density at the center of 
the nucleon is about as high as in the center of Ccl^q. 
Therefore, large number of gluons does not yet imply 
that all of them merge in the transverse plane into more 
or less homogeneous distribution: the positions of the in- 
coming nucleons still remain the dominant source of the 
initial state fluctuations. 

As the "gluonic spots" from single nucleons remain the 
source of initial state fluctuations, one may ask if numer- 



ous partons coming from it may be correlated in their 
quantum numbers, forming specific large- amplitude color 
fields. One particular kind of such fields got special atten- 
tion: those are the topologically nontrivial gluonic con- 
figurations called the QCD sphalerons [25]. They are glu- 
onic field configurations which originate from excitations 
of the topologically nontrivial vacuum fields (instantons) . 
While they rapidly explode into multiple gluon state, 
they strongly violate CP and chiral symmetries locally, 
producing in particularly ±2Nf ( « 6) units of chiral- 
ity per sphaleron. As pointed out in [26 , such strong 
local CP violation induces special event-by-event fluctu- 
ations in the CP-odd observables, e.g. they should induce 
charge asymmetry along the magnetic field. Clearly those 
should be looked at in special studies. 

Another coherent color field configurations which de- 
serve to be specially studied are (colorelectric) flux tubes. 
In pp collisions the view that a field created by longitu- 
dinally separated charges makes a flux tube is a conse- 
quence of confinement, and thus must happen in vac- 
uum. Many popular event generators are based on the 
Lund model, cleverly parameterizing flux tube produc- 
tion and decay. In AA low energy collisions many flux 
tubes are produced, and their possible fluctuations into 
the so called "color ropes" has been studied, initiated by 
the paper of [22]. If two elementary color charges can be 
combined, they may either cancel each other or produce 
higher representations of the SU(3) group, in which case 
the rope energy (and entropy, after its decay) is propor- 
tional to its flux squared. Further applications of these 
ideas for strangeness production in AA collisions can be 
found in [23] . 

Studies of the flux tubes lay dormant till recent dis- 
covery of the so called "hard ridge" [28] by STAR collab- 
oration at RHIC. One possible origin of it [13 is hydro- 
carried longitudinal flux tube, created at the hard colli- 
sion point. This explanation may work provided the flux 
tubes survive as such till freezeout: as was argued in that 
paper this indeed should happen at the periphery of the 
fireball, where matter is not far from the deconflnement 
transition, forming a kind of "dual corona" of the QGP 
fireball similar to Solar corona full of flux tubes. 



IV. WHICH CORRELATIONS SHOULD ONE 
STUDY? 

A. Central collisions: two versus many-particle 
correlators 

Suppose a given event has certain (2-dimensional) dis- 
tribution over transverse momenta of the secondaries. 
This distribution can be decomposed into Taylor series 
of the momenta p x -,Py or into the angular harmonics 



dN 
dptd<t> 



f(Pt 



1 + (2a n cos(n(j)) + 2b n sin(n<j))) 



n>0 



a n =< cos{n(j)) >,6 n =< sin(n(j)) > (4.1) 
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(Note that instead of the square bracket one can also 
use the exponent of the sum, which will include trivial 
higher order effects and enforces positivity of the distri- 
bution: but we would assume all v n to be very small, 
for simplicity.) Instead of using the a, b pair, one may 
also introduce the modula and the phases writing it as 
2v n cos[n((j) — £ n )] with positive v n . In order to simplify 
subsequent formulae, we however prefer to write it using 
the complex exponent 



dN 

dp 2 t d(j) 



fiPt 



(4.2) 



where the sum goes over all integer n, positive and neg- 
ative, with v$ = 1 and V- n = v n . 

The sum of the harmonics presumably describe some 
interesting shapes. For example, those can be the sound 
circles from point perturbations, leading to two-maxima 
shapes mentioned in the Introduction. Unfortunately, 
there are multiple perturbations in any event, and the 
fluctuations we discuss are so small that they can only 
be studied by finding statistically significant small corre- 
lations, using large (~ 10 9 ) sample of available events. 

Discussing how to do it let us, for simplicity, first dis- 
cuss a subset of exactly central collisions, with zero (neg- 
ligibly small) impact parameter and thus exact azimuthal 
symmetry of the event ensemble. Individual events would 
have of course some fluctuations pointing to some direc- 
tions £ n , but their absolute directions are arbitrary, not 
defined relative to anything physical. Whatever fluctua- 
tion happens, it is clear that its copy rotated by any angle 
should also exist and have the same probability. Say, if 
the elementary perturbation is local (delta-function-like 
in the transverse plane): its angular position in the trans- 
verse plane is the only meaningful azimuthal orientation. 
Let us express it as 



£n — £,p H~ <^n 



(4.3) 



where the tilde indicates the angle relative to the pertur- 
bation. Obviously £ p is a random variable, changing from 
event to event, and thus should be averaged out. It turns 
out that the way the correlations work out is quite dif- 
ferent for (i) the two-body and (ii) the many-body (three 
or more) correlation functions. Indeed, in order to get 
the two-body correlation function one has to take a the 
square of the single-body distribution (4.2) 



^ v^v^explin^x + m 2 2 

n± ,ri2 

-iniim ~ in 2in 2 ~ i(ni + n 2 )£ p ] 

and average it over £ p randomly distributed over the cir- 
cle. As a result only terms satisfying 



amplitude e\. Second, it becomes a function of the an- 
gular dif ference A0 = 0i — 02, as expected. And, last 
but not least, all the phases £ n disappear. 

This facts are of course well known , usually expressed 
as harmonics of the 2-body correlator 



C 2 (A0) =< 



d 2 N 



> 



CnA 



/ d{A4>)C 2 {A4>)cos{nA4>) 
Jd(A<f>)C 2 



=< vl > 



(4.5) 



(4.6) 



So, the two-body correlator provides squared amplitudes 
of the original harmonics, averaged over the events. This 
is e.g. how Alver and Roland [12] and others have ob- 
tained their estimates for the "triangular" flow. 

However, the situation is different for manybody (three 
or more) corr elat ion functions. Indeed, if the single-body 
distribution ( |4.2[ ) is cubed (or raised into higher power), 
one finds a triple sum in which random perturbation di- 
rection appears as exp[i(ni-\-n2-\-ns)i/j p ]. Averaging over 
it leads now to the condition 

m + n 2 + n 3 = (4.7) 

Eliminating e.g. ns one finds the double sum 



ri\ ,77,2 



^ri2 ^ni 

+ n2 exp{i[rii((j)i — 03) + n2(02 — 03) 

— ™l(£m - £ni+n 2 ) - ™2(£ni £ni+n 2 )]} 



in which the relative orientations of the different harmon- 
ics are still present. Taking corresponding moments over 
the observed angle differences 0i — 03,02 — 03 one can 
single out the corresponding terms, and thus measure the 
combination 

< V ni V n2 V n3 COs(ni£ ni + n 2 £n 2 + ^3^n 3 ) > ( 4 -8) 

in which t hree integers are subject to the "resonance" 
condition (4.7). Therefore, since there are many mo- 



n\ + n 2 = 



(4.4) 



ments of various many-body correlators, much more com- 
plete data analysis is possible, ultimately allowing to ob- 
serve also the phases of the harmonics. The price to pay 
is that the smallness of all harmonics now appear in the 
third (or higher) power, so it is more difficult to get its 
values out of the statistical noise. 

Let us now briefly discuss how comparison to theory 
should be done, assuming such averages are experimen- 
tally measured. There are two steps to be done. First, 
using linear response hydro one can approximate it as 

< V ni Vn2Vn 3 COs(ni£ ni + n 2 £n 2 + ^3^n 3 ) > ( 4 - 

( — )( — )( — ) < e ni e n2 e n3 cos(ni£ni + U 2^n 2 +^3^n 3 ) > 



survive, which has three important consequences. First, 
a double sum collapses into a single sum with the squared 



Second, one should c hang e from the flow angles to defor- 
mation angles using (3.2). Note that in each three terms 
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Ui in numerator and denominator cancel, leaving only 3tt 
or simply the sign change 

cos(ni£ ni + n 2 £n 2 + n 3 ^ n3 ) = (4.10) 

-COs(niVv + n 2^n 2 + ™3^n 3 ) 

The resulting correlation of the amplitudes and orienta- 
tions of the initial state fluctuations can be calculated 
from initial state model, as we have done above for the 
first three harmonics resonance 1 + 2 = 3 as well as 
3 + 2 = 5. 

In fact it is not necessary to take angular moments of 
the correlation functions. Instead of lengthy but obvious 
general expressions, let us give a simple but instructive 
example. We are interested in telling the difference be- 
tween the two-prong events predicted by sound propaga- 
tion and three-symmetric case of purely triangular flow. 

Consider the two simplest shapes depicted in FigfTTj 
The upper figure has (green) curve which shows the two- 
peak distribution, and the lower figure has (red) three- 
peak one. If it is projected into angular harmonics, and 
the (viscous) filter kills all but the first three, one gets 
their approximation by the other lines. The 1-st and 3-ed 
harmonics are correlated because of absent third peak. 

Defining the 3-particle correlator as for the 2 particle 
one, with the averaging over the perturbation angle 

d 2 N 

C3( ^'^^ )=< #I#^ > ^ (411) 

One finds it to be a function of two angles, say <pi — 
</>3, 02 — </>3- The corresponding plots for two examples in 
question, the 2-peak and 3-peak distributions, is shown 
in Fig(T2] Comparing two plots one can see how different 
they are. In experiment, which will produce a mixture 
of the two, one would see that intermediate pattern, and 
hopefully determine the presence and probability of both 
type of events. 

B. Mid- central collisions and the two-body 
correlators relative to the event plane 




FIG. 11: (upper) Two-peaked (green) curve shows the sim- 
plest model of angular distribution to be discussed. The (red) 
oscillating curve is its version after the filter removes all but 
harmonics with n > 3. 

(lower) The lower (red) curve shows the dependence of the 2- 
particle correlation function on cj)\ — 02 for the filtered 2-peak 
distribution. The higher (green) curve is shown for compari- 
son, it corresponds to another model distribution with three 
identical peak shifted by 120°. 



Nonzero impact parameter violates axial symmetry 
and creates "directed flows" , e,g, the famous elliptic flow 
with nonzero < v 2 > z / z 0. By mid-central collisions we 
mean a centrality region in which < e 2 > is large itself, 
and is also large compared to its fluctuations (recall that 
it is not so for central and very peripheral collisions). For 
example, e 2 is 0.5 (0.3) for N p = 100 (200) participants, 
with Se 2 ~ 0.1. Furthermore, as seen in Figj6}}, its angle 
ijj 2 is very much directed a t + 7r/2 (the tips of the al- 
mond) and therefore (using ( |3.2| for n = 2) the flow angle 
is peaked "in-plane" ,£2 = 0, 7r, as indeed observed. 

If so, for one of the harmonics being the second e.g. 
ns — ±2 one can approximate a product of three defor- 
mations as follows 

< v nx v n ^v n?) cos(n x i nx + n 2 £n 2 + n 3 £ n3 ) > (4.12) 
U71 ^2 

_ (^l)(^-)(_ ) < e 2 >< e ni e n2 co,s(ni^ ni +n 2 £ n2 ) > 



by factoring out large and non-fluctuating < 62 > from 
two other harmonics which are small and fluctuating. 
Note that the resonance condition now means ri\ zb 2 = 
n 2l and that by putting £2 = we have selected a frame 
in which the (experimentally determined reaction plane) 
is the x axis. 

Basically the lesson here is that for mid-central col- 
lisions the "reaction plane" plays the same role as the 
third body, so we are reduced to two small fluctuating 
and correlated harmonics. The simplest nontrivial ex- 
ample of resonance condition of the kind is 3-1-2=0 (re- 
cently studied by Teaney and Yan [27]). while the next 
is 5-3-2=0. 

We had already calculated the combinations of two 
fluctuating harmonics with the appropriate cosines 
above, for these two cases, in the Glauber model. They 
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FIG. 12: The contour plot of the three particle correlators 
as a function of <pi — 03 and 02 — 03, for 2-peak distribution 
(upper figure) and 3-peak one (lower). As one can see, their 
shapes are quite different. 



are by no means small: for example for the centrality 
class with 100 participants they are 

< eie 3 cos (3^3 - ^i) >~ -0.015 (4.13) 

< e 3 e 5 cos (3^3 - 5^ 5 ) >~ -0.05 (4.14) 

and therefore we expect it to be observable, with about 
as large statistics as needed for the usual quadratic fluc- 
tuations. 



V. SUMMARY 

In this work we have (i) discussed the setting, iden- 
tifying the main scales of the problem. Then (ii) we 
studied in detail the initial state fluctuations originat- 
ing from nucleon positions, emphasizing existence of the 
nontrivial phase relations between different harmonics. 
Finally, (iii) we discussed correlations functions of 2 and 
many hadrons, pointed out the principle difference be- 
tween them, with the latter allowing to experimentally 
measure the relative phases of these harmonics. 



Let us now recapitulate the lessons from this study in 
a bit more detail. 

Unfortunately, the perturbations we speak about are 
too small to be measured directly, on event-by-event ba- 
sis, and should be instead reconstructed from the statis- 
tically obtained correlation functions. One good thing 
coming from it is that multiple but independent fluc- 
tuations from local perturbations in a single event and 
in the ensemble are treated in one and the same corre- 
lation function, while all trivial effects are statistically 
subtracted. Traditionally the initial state perturbations 
and final state corrections to collective flow are consid- 
ered in a form of their angular harmonics, which we call 
e n and v n respectively. Their ratios v n /e n can be cal- 
culated by the linearized hydrodynamics: the details of 
that is the subject of our next paper. 

"Harmonic flows" have been so far treated as indepen- 
dent or incoherent fluctuations, added in quadratures. 
We however pointed out that it is only true for say cen- 
tral collisions and two-particle correlations. For central 
collisions many (3 and more) particle correlations already 
include phases of the deformations i/j n , as do two-body 
correlations relative to event plane for mid-central colli- 
sions. Many different moments related to "resonances" 
between 3 integers can be measured, providing experi- 
mental opportunity to find out how the perturbations are 
correlated and eventually to refine models of the initial 
state. The magnitude of such terms, as we have shown for 
Glauber model, is no smaller than for the terms already 
studied. 

Coherence in phases of the deformations imply the 
interferences between the harmonics of the flow. Only 
adding them together one can follow how small "hot" (or 
"cold") spots created by quantum fluctuations of inter- 
acting nucleons propagate away from the point of origin. 
Only in this way one can understand the role played by 
the "hydrodynamical causality" , insisting that large part 
of the fireball should remain completely unaffected by 
the perturbation since the signal cannot possibly reach 
it before the freezeout. Only a complete Green func- 
tion, collecting all hydro harmonics, will describe correct 
shapes of propagating waves, such as e.g. (i) the "cylin- 
ders" from the initial state fluctuations, or (ii) the "Mach 
cones" from quenched jets. 

There are two basic scales defining those perturbations, 
the sound hori zon H s (2.2) and the "the viscous horizon 
scale" R v (2.7). The former gives the size of the pertur- 



bation, stemming from a local perturbation, the second 
its width. We have for example argued that by chang- 
ing the centrality of the collisions, one can change the 
relation between the (smaller) fireball size and the sound 
horizon: this should dramatically change the shape of the 
event (see Fig. 4 for explanation). It is an important ob- 
jective of the experimental heavy ion program in general 
to measure these two scales, extracting experimental val- 
ues of the speed of sound and viscosity. It can be done for 
example by changing centrality and observing the change 
of shape of the underlying event. 
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